Maternal carryover, winter severity, and brown bear abundance relate to elk demographics

Ungulates are key components of ecosystems due to their effects on lower trophic levels, role as prey, and value for recreational and subsistence harvests. Understanding factors that drive ungulate population dynamics can inform protection of important habitat and successful management of populations. To ascertain correlates of ungulate population dynamics, we evaluated the effects of five non-exclusive hypotheses on ungulate abundance and recruitment: winter severity, spring nutritional limitation (spring bottleneck), summer-autumn maternal condition carryover, predation, and timber harvest. We used weather, reconstructed brown bear (Ursus arctos) abundance, and timber harvest data to estimate support for these hypotheses on early calf recruitment (calves per 100 adult females in July–August) and population counts of Roosevelt elk (Cervus canadensis roosevelti) on Afognak and Raspberry islands, Alaska, USA, 1958–2020. Increasing winter temperatures positively affected elk abundance, supporting the winter severity hypothesis, while a later first fall freeze had a positive effect on elk recruitment, supporting the maternal carry-over hypothesis. Increased brown bear abundance was negatively associated with elk recruitment, supporting the predation hypothesis. Recruitment was unaffected by spring climate conditions or timber harvest. Severe winter weather likely increased elk energy deficits, reducing elk survival and subsequent abundance in the following year. Colder and shorter falls likely reduced late-season forage, resulting in poor maternal condition which limited elk recruitment more than winter severity or late-winter nutritional bottlenecks. Our results additionally demonstrated potential negative effects of brown bears on elk recruitment. The apparent long-term decline in elk recruitment did not result in a decline of abundance, which suggests that less severe winters may increase elk survival and counteract the potential effects of predation on elk abundance.


Introduction
Ungulates are important components of ecosystems due to their effects on lower trophic levels and role as prey, and also have economic and cultural value because of human recreational and subsistence harvests [1,2]. Ungulate abundance and recruitment estimates are commonly used to inform population (e.g. harvest quotas) and land (e.g. timber harvest scheduling) management decisions [3,4]. Understanding factors that drive ungulate abundance and recruitment is similarly important, as it further informs managers' ability to protect important habitat and manage populations [5].
Numerous factors can influence ungulate population dynamics including weather, predation, forage availability and quality, and interactions among these factors [6,7]. Determining which factors drive ungulate population dynamics is especially important with current unprecedented changes in species' ranges and population sizes due to anthropogenic habitat alteration and climate change [8,9], with wide-reaching implications for ecosystem and wildlife management [10]. Long-term data provide opportunities to draw stronger inference which can improve our understanding of factors that influence ungulate recruitment and population dynamics [3,8].
Weather can affect ungulate population dynamics directly through temperature-driven physiological stress, and indirectly by altering forage availability [9]. Winter severity can reduce maternal ungulate body condition and juvenile survival through increased thermoregulatory costs from cold temperatures [11], and increased difficulty of movement or reduced forage availability due to increased snowfall [12,13]. Spring weather also can affect ungulate population dynamics [14] because later, colder, and drier springs limit vegetation growth and availability [15]. Reduced forage quantity and quality can cause a spring nutritional "bottleneck" that may cause ungulate energetic costs to exceed reserves, reducing survival and reproductive success [16,17]. Similarly, summers or falls with low temperatures and precipitation can reduce summer and fall forage availability, in turn reducing ungulate nutritional carryover, which reduces pregnancy rates, winter survival, body condition, and ultimately abundance and recruitment [18,19]. Predators also affect ungulate survival [20,21], which influences long-term abundance and recruitment [3,22]. Anthropogenic habitat alterations such as timber harvest influence ungulate habitat suitability, and therefore abundance and recruitment [23,24]. Recently harvested timber stands can increase habitat suitability by increasing forage availability [4,25] and young stands can increase habitat suitability by increasing hiding cover [4], however timber harvest also reduces mature forest area, which may provide shelter from predation and severe winter weather [26,27].
To quantify factors that limit ungulate populations, we tested five non-exclusive hypotheses to assess their effects on Roosevelt elk (Cervus canadensis roosevelti) early calf recruitment (defined as calves per 100 adult females in July-August; henceforth recruitment) and abundance ( Table 1). The winter severity hypothesis predicts that cold temperatures or high snowfall cause decreased abundance and recruitment [12,21]. The spring bottleneck hypothesis predicts that a prolonged winter or a spring that is cool or dry leads to decreased abundance and recruitment [17,28]. The maternal carryover hypothesis predicts that a short, cool, or dry summer; a cool or dry autumn; or an early first freeze causes decreased abundance and recruitment the following summer [18]. The predation hypothesis predicts that increasing abundance of brown bears will decrease elk recruitment, as brown bears may predate elk calves [3,19,21]. Finally, the timber harvest hypothesis predicts that timber harvest will alter elk habitat suitability, and therefore recruitment [4,29], in that an increase in area of recently harvested timber stands will increase forage availability and therefore positively influence recruitment, or alternatively, an increase in area of recently harvested timber stands will reduce availability of mature forest and therefore negatively affect recruitment [4].

Study area
Afognak (1,809 km 2 ; 58.3279˚N, 152.6415˚W) and Raspberry (197 km 2 ; 58.0708˚N, 153.1876˚W) islands are in the Kodiak Archipelago, Alaska, USA (Fig 1). The islands are 5 km north of Kodiak Island and separated by a 1.5-km wide strait. Afognak and Raspberry islands are primarily owned by Native corporations (64%), followed by state (27%) and federal (9%) ownership. Both islands contain gradual sloping mountains ranging from 300 to 800 m in

Elk data
To measure elk early calf recruitment and abundance, we used aerial surveys conducted sporadically during 1958-1964 and annually thereafter [35]. Elk counts were conducted by ADFG biologists twice per year, first during July-August, emphasizing on herd composition to measure early calf recruitment (calves per 100 adult females; age ratios), and then during September, to estimate population size [35]. Surveys were conducted with no set flight plan at various altitudes to maximize elk sightability and identification, and counts were conducted independently by the pilot and biologist until a consensus was reached [35]. For composition counts, adult females, males and calves were identified by size, antler presence, and coloration [44]. During 1992-2020, locations of elk herds for surveys were obtained with very high frequency radio-or Global Positioning System (GPS) collars deployed as part of a parallel study [35].
While extensive canopy cover limited efficacy of herd abundance counts on Afognak Island, composition counts only require a portion of the herd to be observed, with the assumption that the portion observed is representative of the herd, and therefore can be conducted in areas where population counts are not feasible [45]. Because Raspberry Island has little canopy cover, we were able to incorporate elk abundance for the Raspberry herd only. We accounted for elk harvest on Raspberry Island when modeling elk abundance. Elk harvest occurs during fall after population counts (e.g. 25 September-30 November in 2021) and harvest information has been collected by ADFG since 1949 [42].

Climate data
We used weather data from Kodiak Airport [30,46], 35 km south of Afognak Island, to derive predictor variables to test for support of the winter severity, spring bottleneck, and maternal carryover hypotheses (Table 1). For the winter severity hypothesis, we tested mean temperature, mean monthly minimum temperature, total snowfall, and mean daily snow depth during winter (1 November-31 March). For the spring bottleneck hypothesis, we tested mean spring (1 April-31 May) temperature, total spring precipitation, date of last spring freeze, and the number of days between the first fall freeze and last spring freeze (i.e., winter duration). For the maternal carryover hypothesis, all predictors were from the summer and autumn in the year before the respective population count and composition surveys. We tested mean summer (1 June-31 August) temperature, total summer precipitation, growing degree days calculated at base 5˚C [47], the standardized precipitation evapotranspiration index (SPEI) [48], mean fall (1 September-31 October) temperature, total fall precipitation, and date of first fall freeze.
In boreal regions, SPEI is negatively correlated with summer productivity [49]. We used R package SPEI to calculate potential evapotranspiration using the "thornthwaite" function [50] and "spei" to calculate SPEI with a calibration period during January 1950-May 2021 [48]. We used two SPEI calculations: a 6-month period (January-July) which best predicts forest productivity (forest SPEI), and a 3-month period (July-September) which best predicts grassland productivity (grassland SPEI) [51].

Brown bear data
Brown bear harvests on the Kodiak Archipelago occur during fall and spring (25 October-30 November and 1 April-15 May) [42]. Hunters are limited to one bear per permit and harvest of females with dependent young is prohibited. All harvested bears are inspected by ADFG personnel who record sex and location of kill for each bear and extract a vestigial premolar tooth for cementum aging [52]. Additionally, ADFG inspects and extracts teeth for all other known brown bear mortalities (e.g. illegal harvest, defense of life and property, vehicle or natural mortalities).
We used Downing population reconstruction to estimate brown bear abundance on Afognak and Raspberry islands using total harvest by year and harvest-by-age data [53,54]. Downing population reconstruction assumes that the primary source of mortality is harvest, but we augmented our data by including all known mortalities [53]. Downing reconstruction also implies that cohort harvest mortality relative to total mortality is constant among years, that mortality rates for the oldest two reconstructed age classes are equal, and that the aged sample is unbiased, but Downing's method is robust to violations of many of these assumptions [54]. Additionally, because this brown bear population has a low rate of natural mortality with little interannual variability [41], a high harvest reporting rate [55] and a high and unbiased rate of aging of harvested bears (average number of teeth aged relative to bears harvested = 0.92), we believe that the Downing population reconstruction provided reasonable estimates of brown bear abundance trends.
Cementum aging records were available from 1967 to 2017 (S1 Dataset). Though brown bears do not achieve full body size until 8-14 years old [56], Downing population reconstruction is robust to collapsing age classes [54]. Downing population estimates calculated from these data were similar using 5, 7, and 12 age classes. We therefore used five age classes to increase the number of years available for analyses. As Downing population reconstruction provides inaccurate abundance estimates for one less than the number of age classes estimated in most recent years [53], we excluded the four most recent years of abundance estimates (2014-2017) from our analysis. We did not analyze sexes separately because a preliminary Downing reconstruction using all known mortalities of bears in the Kodiak Archipelago revealed no difference in abundance estimates when sexes were analyzed separately (average yearly proportional difference in summed sex-separated estimate compared to sexes combined = 0.022).

Timber harvest data
We compiled 45 years (1976-2020) of timber harvest data from Afognak Native Corporation, Koniag Native Corporation, Koncor, Natives of Kodiak Native Corporation, and Ouzinkie Native Corporation. We verified, corrected, and added to these data using Google Earth and Landsat satellite imagery [57,58]. We characterized forest stand by age based on their habitat value for elk including clear-cut (<1 year), early regeneration (1-5 years), late regeneration (6-30 years), and mature forest (>30 years) [59][60][61], then calculated the total area harvested and area of each age class in each year (Table 1).

Analysis
Using location data from GPS-collared elk [39], we separated recruitment data based on elk herds which occurred in areas with timber harvest (Marka Lake, Waterfall, Portage Lake, Duck Mountain, Seal Bay; "timber harvest dataset") and areas without timber harvest (Raspberry, Malina; "non-timber harvest dataset"; Fig 2) by summing calf and cow counts from the respective herds then calculating the number of calves per 100 adult females (hereafter age ratio). We also used a similarly calculated combined dataset of all herds except Tonki as dense forest cover precluded accurate composition counts ("island-wide dataset"). Before summation of the recruitment datasets, we removed herd counts for years in which the number of females was less than the number of calves, because this indicated an error in data recording or incomplete count of the grouping; the number of calves counted was zero, because these indicated that the survey was conducted too late in the year to distinguish cows from calves; and years in which cows and bulls were not counted separately (S1 Table). We excluded the recruitment count in 2006 from the timber harvest dataset because the total number of elk counted across all herds was five. Given these exclusions, we did not analyze recruitment data for 1998, 2006, or 2007 (S2 Dataset). We standardized all predictor variables by mean and standard deviation to compare effect sizes [62].
We fit linear models to each elk dataset (timber recruitment, non-timber recruitment, island-wide recruitment, and Raspberry herd abundance) with year as the predictor, and used p-values (α < 0.05) to determine long-term trends [62]. To determine long-term trends in estimated brown bear abundance, we also fit linear models to the brown bear population reconstruction and natural mortality rates, also using p-values (α < 0.05).
The winter severity, spring bottleneck, maternal carryover, and timber harvest hypotheses were represented by several intercorrelated predictor variables that had similar biological meanings (Table 1; S1 File). To determine which measured predictor variable best represented the corresponding hypothesis, we first fit univariate regression models for each predictor to determine which had the strongest relationships with the response variables (i.e. recruitment and abundance). Within each hypothesis, we then selected the predictor with the greatest R 2 value to represent the respective hypothesis in the second stage of the analyses [62,63]. We examined the correlation between selected predictors using Pearson's product-moment correlation coefficient (r). For those highly correlated (|r| > 0.7), we replaced the correlated predictor with an uncorrelated predictor with the next greatest R 2 value that allowed us to test the maximum number of hypotheses [62,63].
To test our hypotheses on recruitment, we used the selected predictors for each hypothesis to fit Bayesian general linear models to the timber harvest recruitment, non-timber harvest recruitment, and island-wide recruitment datasets (S1 Appendix). Area of timber harvest was included only in the timber harvest and island-wide recruitment models.
To test our hypotheses on changes in Raspberry herd abundance, we used the selected predictors for each hypothesis to fit a Bayesian Gompertz state-space model to the time series of elk count data, which accounts for potential density dependence and has been used previously to assess drivers of ungulate population dynamics [7,12,64]. We also explicitly incorporated known harvest size in each year [64]. We did not include reconstructed brown bear abundance as a covariate for this model because we were unable to accurately estimate brown bear abundance on Raspberry Island. Population count data were imputed using linear interpolation for 1959, 1960, 1963. A more detailed description of the model can be found in S1 Appendix. Because the coefficient for density dependence in the Gompertz model with external predictors is only partially identifiable even within a state-space context [65,66], we also tested for density dependent effects within the Raspberry herd using the Dennis-Taper parametric bootstrap likelihood ratio t-test, and examined a plot of the natural log of the per capita growth rate versus abundance [66,67].
We used Bayesian methods for all models because they allow for clear interpretation of results and have better inference with smaller sample sizes [68]. We fit models using Gibbs sampling through a wrapper for the R package rjags, jagsUI (annotated jagsUI code provided in S2 File) [69,70]. We used uninformative normally distributed (for parameters) or uniformly distributed (for error) priors with three chains, a thin rate of 5, and 50,000 iterations with a burn-in of 10,000. We examined traceplots, residual plots, scale reduction factors (R), and Bayesian p-values to assess model convergence and goodness of fit, respectively, and used 95% credible intervals to determine influence of predictors [68].
We used 39 years of age ratio data during 1967-2013 from the Waterfall, Duck Mountain, Marka Lake, Portage Lake, and Seal Bay herds to assess effects of timber harvest on elk recruitment. The elk age ratio ranged from 16 to 48 calves per 100 adult females, with a mean age ratio of 33±7.9 calves per 100 adult females and demonstrated no long-term trend (estimate = -0.15, R 2 = 0.043, p = 0.104, t = -1.664; S2 Fig). For the timber harvest hypothesis, the two predictors with the greatest R 2 values were correlated with brown bear abundance (timber harvest aged 6-30, r = 0.955; total area harvested, r = 0.927; Table 2; S1 File), therefore we used area of timber harvest aged 1-5. None of the other selected predictors (i.e. those with the greatest R 2 values for each hypothesis) were highly correlated (all |r| < 0.7) and therefore were used in the final model. The final model included mean winter temperature representing the winter severity hypothesis, winter duration representing the spring bottleneck hypothesis, day of first fall freeze representing the maternal carryover hypothesis, brown bear abundance representing the predation hypothesis, and area of timber harvest aged 1-5 for the timber harvest hypothesis. The final model successfully converged (allR < 1.1) and fit well (Bayesian p- PLOS ONE value = 0.538), but there was no support for any of the predictor variables, as all 95% credible intervals overlapped zero (Table 3). We used 44 years of age ratio data during 1967-2013 from the Malina and Raspberry herds to assess drivers of elk recruitment in areas without timber harvest. The elk age ratio ranged from 16 to 57 calves per 100 adult females, with a mean ratio of 33±9.4 calves per 100 adult females (S2 Fig). The ratio decreased by 0.32 calves per 100 adult females annually (R 2 = 0.185, p = 0.001, t = -3.316). None of the selected predictors (i.e. those with the greatest R 2 values for each hypothesis) were correlated (all |r| < 0.7; S1 File) and therefore all were used in the final model. The final model included mean monthly minimum winter temperature representing the winter severity hypothesis, mean spring temperature representing the spring bottleneck hypothesis, day of first fall freeze representing the maternal carry-over hypothesis, and brown bear abundance representing the predation hypothesis ( Table 2). The final model successfully converged (allR < 1.1) and fit well (Bayesian p-value = 0.528). Brown bear abundance was the only influential predictor, as all other predictors had 95% credible intervals overlapping zero  Fig 3). We used 44 years of age ratio data during 1967-2013 from all elk herds in the island-wide recruitment model. The elk age ratio ranged from 20 to 57 calves per 100 adult females, with a mean ratio of 33±7.2 calves per 100 adult females. The ratio decreased by 0.22 calves per 100 adult females annually (R 2 = 0.165, p = 0.003, t = -3.116; S2 Fig). The two predictors with the greatest R 2 values for the timber harvest hypothesis were correlated with brown bear

PLOS ONE
abundance (timber harvest aged 6-30 years old r = 0.959, total area harvested r = 0.929; Table 2; S1 File), so we used area of timber aged 1-5 for the timber harvest hypothesis. None of the other selected predictors (i.e. those with the greatest R 2 values for each hypothesis) were correlated (all |r| < 0.7) and thus were used in the final model. The final island-wide recruitment model included mean monthly minimum winter temperature representing the winter severity hypothesis, winter duration representing the spring bottleneck hypothesis, day of first fall freeze representing the maternal carry-over hypothesis, brown bear abundance representing the predation hypothesis, and area of timber harvest aged 1-5 for the timber harvest hypothesis. The final model successfully converged (allR < 1.1) and fit well (Bayesian pvalue = 0.532). Brown bear abundance and date of first fall freeze were the only influential predictors, as other 95% credible intervals overlapped zero ( Table 3). The predicted elk age ratio increased as first fall freeze occurred later in the year, increasing by 14.3 calves per 100 adult females across the range of observed values (17 September-27 October; Fig 4A). As brown bear abundance increased, the predicted elk age ratio decreased by 11.2 calves per 100 adult females across the range of observed values (83-331 individuals; Fig 4B). The negative effects of increasing brown bear abundance on elk recruitment (β = -3.42, 95% credible interval = -5.89, -1.07) were greater than the positive effects of a later fall freeze (β = 2.80, 95% credible interval = 0.21, 5.47).
We used 63 years of aerial count data during 1958-2020 for the Raspberry herd state-space abundance model. Raspberry herd population counts ranged from 10 to 230 elk with a mean of 116±61.7 elk (mean density = 0.59±0.31 elk/km 2 ), and abundance appeared cyclical with no evidence of long-term trend (estimate = 0.723, R 2 = 0.026, p = 0.117, t = 1.594; S3 Fig). None of the selected predictors (i.e. those with the greatest R 2 values for each hypothesis) were correlated (all |r| < 0.7; S1 File) and therefore all were used in the final model. The final model for growth rate included mean winter temperature representing the winter severity hypothesis,

PLOS ONE
mean spring temperature representing the spring bottleneck hypothesis, mean fall temperature representing the maternal carryover hypothesis, and a parameter for density dependence ( Table 2). The final model successfully converged (allR < 0.01) and fit well (Bayesian pvalue = 0.528). Mean winter temperature and the parameter associated with density dependence were both influential predictors; all others had 95% credible intervals overlapping zero ( Table 3). As mean winter temperature increased, predicted growth rate increased by 0.53 over the range of observed values (-2.80-2.93˚C; Fig 5A). Though the credible interval of the density dependence parameter did not overlap zero (Fig 5B), the Dennis-Taper parametric bootstrap likelihood ratio t-test was insignificant, indicating lack of density dependence (p = 0.336; S4 Fig).

Discussion
Weather explained variation in elk abundance and recruitment. The winter severity hypothesis was supported by the Raspberry herd abundance model, while the maternal carryover hypothesis was supported by the island-wide recruitment model. We found no support for the spring bottleneck hypothesis, but we found the predation hypothesis received the greatest overall support in explaining recruitment island-wide and in areas without timber harvest. We found no support for the timber harvest hypothesis.
We found that elk abundance was affected by winter temperature, supporting the winter severity hypothesis. Increased winter severity, in the form of increased snowfall [21], snow depth [13], and decreased temperatures [71] can have a negative effect on elk survival across a broad range of climatic conditions, which in turn leads to decreased abundance after harsh winters [72]. Though severe winters and later, colder springs can be detrimental to ungulate recruitment [14,73,74], our findings on early calf recruitment did not support the winter severity or spring bottleneck hypotheses, indicating that winter severity has an impact on overwinter survival, rather than successful parturition and survival of calves to two months.. Winter may have less impact on early calf parturition for northern ungulates in coastal areas with milder winters and warmer summers than typical mainland subpolar and continental climates [6,75]. For example, caribou early calf recruitment was best predicted by higher availability and quality of forage on summer range than winter ranges on the Kenai Peninsula, Alaska [76], indicating support for the maternal carryover hypothesis rather than the spring bottleneck or winter severity hypotheses.
We found that a later fall freeze positively impacted island-wide recruitment, supporting the maternal carryover hypothesis. The positive effects of higher quality forage in summer and fall due to increased precipitation and warmer temperatures on ungulate abundance and recruitment are well-documented [19,72,76]. However, little research has been conducted on how winter onset timing affects ungulate recruitment [77]. Though fall forage quality decline can be gradual, it is accelerated by sudden frost events [78]. Increased availability of fall forage over a longer time may be crucial for ungulates [79], as high-quality fall forage can increase ungulate body fat levels at the beginning of winter, therefore increasing winter survival, successful parturition, and spring calf condition [80,81]. Additionally, warmer or longer falls allow later born calves to increase their body condition to catch up with earlier-born conspecifics, which in turn can increase their winter survival and likelihood of pregnancy as yearlings [18].
Elk recruitment declined with increasing brown bear abundance in the non-timber harvest and island-wide recruitment models, supporting the predation hypothesis. Predation by brown bears on ungulates, specifically neonates, is well documented, as are the effects of brown bear presence on elk reproductive success and calf survival [3,20], but few studies have examined the effects of brown bear abundance on ungulate recruitment, especially in a northern system devoid of other predators. Though brown bears can be specialist predators of ungulate neonates in the spring, predation by brown bears on ungulates is likely limiting rather than regulatory, especially in higher-density ungulate populations [40]. Because brown bears often predate ungulates opportunistically [82], their role as predators varies based on the availability of ungulates and presence of alternate food [40,83]. Brown bear predation patterns are strongly seasonal, especially in coastal populations, as brown bears switch from ungulate calves in the spring to more readily available salmon and berries during summer [84]. Though winter severity and predation may have additive effects on ungulate survival [21], we did not find that the long-term decline in elk recruitment translated to a decline in abundance, indicating that high juvenile and adult survival in mild winters may be supporting compensatory brown bear predation, even in a low-density population [7,72,85,86]. The strong influence of brown bear abundance on elk recruitment indicates that ungulate populations that are small or declining may experience accelerated declines where brown bear populations are increasing [87], especially where other predator species occur, or where winter conditions are more severe [21,86].
Area of timber harvest did not influence elk recruitment, potentially because the positive effects of forage quantity in recent timber harvest areas negated any adverse effects of reduced cover from severe weather and predation [4,25]. Alternatively, because northern Afognak Island, where timber harvest occurs, is more heavily forested than southern Afognak Island and Raspberry Island, the lack of decline in recruitment in the timber harvest dataset and lack of a predation effect identified from the timber harvest model could indicate that mature and old-growth forests are important for elk during parturition to mitigate brown bear predation on neonates [26,88]. During parturition, ungulates select habitat with high cover and high visibility, such as older forest stands [32,88], which may reduce predation risk [89] and therefore decrease the effect of brown bear predation on recruitment. We recommend further research on the impacts of brown bear predation on elk calf survival, especially in areas with timber harvest.
Our conclusions are subject to consideration of some limitations. First, though the Downing population reconstruction accurately estimates population trends, it generally underestimates abundance by approximately 15% [54]. Therefore, our estimates of brown bear abundance are used only as an index of abundance. Downing population reconstruction assumes there is no variability or long-term trends in harvest rates. Though true rates of harvest relative to brown bear abundance are unknown, number of brown bears harvested varied across years and increased overall, which may have reduced precision and caused overestimation of rate of change in abundance [54]. However, the variability and increase in harvest is likely not due to changes in harvest regulations, number of hunters, or reporting rates [55]. Additionally, the estimated increase in brown bear abundance between 1993 and 2008 is consistent with documented brown bear abundance increases on Kodiak Island during that time [41].
Second, though our models fit well, we were unable to evaluate the effects of density dependence on recruitment because accurate elk abundance estimates were unavailable except for the Raspberry herd. Though in the Raspberry herd abundance model, increasing elk abundance appeared to significantly decrease growth rate and therefore abundance the following year, the effects of density dependence may be overestimated due to independent sampling error and thus weak parameter identification [66,90,91]. Additional analysis found no evidence of density dependence for the Raspberry herd, further limiting our conclusions regarding density dependence in this population. The effects of density dependence tend to be minimal in northern ungulate populations when predators are present [7], indicating that the parameter for density dependence may have been overestimated in our model for Raspberry herd abundance. Though Ricker and θ-logistic models may better approximate the effects of density dependence on abundance for species with slower life histories, these methods also require informative priors to provide accurate parameter estimates for density dependence [65,92]. Furthermore, there is little difference in model fit or parameter estimates between Ricker and Gompertz models, especially when a population is below carrying capacity and the species exhibits undercompensatory population dynamics [12,93,94], which were both likely true for elk in our study [35,40,64].
Third, because over-winter and spring survival of yearlings contributes to recruitment, and yearlings are less likely to produce calves than adult females, two-month recruitment is a measure of calf survival as well as yearling and adult fecundity [44,95]. Furthermore, our results on recruitment trends do not necessarily indicate overall abundance trends, as twomonth recruitment does not measure yearling or adult survival which may drive variation in population growth among ungulate populations [96]. Yet, age ratios can be correlated with annual population growth rates, indicating that they do have value as indicators of population dynamics [3,44]. Age class surveys are subject to observation errors due to sightability and misclassification [97] which could have led to bias in estimated relationships. However, we used composition counts conducted by experienced biologists from the same time of year and removed counts that were incomplete or where detection was compromised, reducing these errors [44]. Finally, because the number of harvest permits issued was directly determined by elk population count in the previous year, and number of elk harvested is determined by number of permits issued, elk harvest and elk abundance were intercorrelated. Therefore, though we believe that harvest likely has a significant impact on elk abundance [21,98], we were unable to determine the effects of harvest on elk abundance or rate of elk population growth.
Forage conditions, predation, and habitat are key determinates of ungulate abundance and recruitment [6,7]. Though drivers of ungulate population variation are species-and area-specific [96], that elk abundance increased with increasing winter mean temperatures supports the hypothesis that elk survival, and therefore abundance, is negatively affected by harsh winter conditions [3,21]. Similarly, that recruitment increased with later first fall freeze adds to a growing body of work which indicates summer and fall forage conditions are more important than winter and spring conditions in driving ungulate vital rates [99][100][101][102], especially in subpolar maritime climates [6,103]. However, the effects of interactions between summer-autumn nutrition, winter severity, and predation on ungulate abundance and recruitment remain unclear. Our results demonstrate a balance between the positive effects of sufficient summerautumn nutrition and favorable winter conditions with the negative effects of predator abundance in influencing ungulate population dynamics. If ungulate abundance and recruitment are declining or low, an increase in predator abundance may cause an increase in additive predation mortality if summer-autumn nutritional conditions are poor or winters are severe [72], especially if there is limited alternative food for predators [82]. Summer and fall forage productivity is important for ungulate populations, especially in areas with increasing predator abundance, as high productivity may balance the negative effects of predation on ungulate recruitment and survival [72,104]. Understanding the underlying mechanisms behind population trends enables us to better develop management and conservation strategies that can plan for variable responses of animals to ecosystem change including habitat alterations and predator abundance.  Table. Elk recruitment data exclusion by herd and exclusion reason. Herd indicates the name of the herd surveyed (locations shown in Fig 2), with the number of years where counts of herds were excluded from analysis because no calves were counted (zero calves), the number of calves was greater than the number of cows (calves > cows), and cows and bulls were not